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Using the results of large scale numerical simulations we study the probability 
distribution of the pseudo critical temperature for the three dimensional Edwards- 
Anderson Ising spin glass and for the fully connected Sherrington-Kirkpatrick model. 
We find that the behaviour of our data is nicely described by straightforward finite- 
size scaling relations. 
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I. INTRODUCTION 

A proper phase transition takes place only in the idealised limit of an infinite number of 
interacting degrees of freedom. Although this limit is never realised in the laboratory (let 
alone in numerical simulations), everyday experience suggests that macroscopic samples are 
infinite for all practical purposes. Spin glasses are an exception. The problem lies 

in their sluggish dynamics at the critical temperature and below. The system remains for 
very long times, or forever, out of equilibrium. In fact, letting the system relax for about 
one hour, the spatial size of the glassy magnetic domains is (at most) of the order of one 



a. 



hundred lattice spacings 

It has become clear lately that, in order to interpret experimental data in spin glasses, the 
relevant equilibrium properties are those of systems of size similar to that of the experi- 



mentally achievable coherence length 



, |5[. Phase transitions on finite systems are actually 



crossover phenomena describable through the well known theory of finite size scaling (see 



e.g. However, a conspicuous feature of disordered systems (and most notably, of spin 
glasses) is to undergo strong sample-to-sample fluctuations in many thermodynamic prop- 
erties. It is thus natural to ask questions about the probability distribution, induced by 
the disorder, of the various physical quantities. Typically the size of these fluctuations de- 
creases when enlarging the size of the equilibrated system; if we wish to have hints about 
their possible relevance in experimental systems, it is important to know the rate at which 
fluctuations decrease with system size. This is particularly important if we are to study dy- 

n n 

namical heterogeneities [7j] in spin glasses [5]. In particular, a relevant but elusive physical 
quantity (potentially relevant to analyse dynamical effects close to the phase transition) is 
the finite-system pseudo-critical temperature. Our scope here is to characterise its statistical 
properties in spin glasses. 

This problem has been extensively studied and is well understood for finite-size weakly 
bond-disordered spin models, below the upper critical dimension d up . For a system of size 
N = L d and a disorder sample J, one can define a pseudo-critical temperature T C J (L) as the 
location of the maximum of a relevant susceptibility: this definition is clearly non unique, 
but all sensible definitions lead to the same scaling behaviour as iV — > oo. According to the 



Harris criterion [8[, a major role is played here by the value of the thermal critical 

exponent of the pure system, vp. If vp > 2/d the disorder is irrelevant, the value of v is not 
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modified by the disorder (i.e. v = z/p), and the width AT C (L) of the probability distribution 
of the pseudo-critical temperature, defined as AT C (L) 2 = E(T^(L) 2 ) — E(T^ (L)) 2 , where 
E(- ■ ■) denotes the disorder average, behaves as AT C (L) oc \jL d l 2 as expected naively 14]. 
In such a situation in the infinite volume limit the (disorder induced) fluctuations of T C J (L) 
are negligible with respect to the width of the critical region and to the finite-size shift of 

T r , that both behave like L^ x l Up . In the other case, when vp < 2 Id. disorder is relevant, the 

J I 

value of v for the disordered model is different from vp and obeys [15( the bound v > 2/d. 
In this case AT C (L), the width of the critical region, and the finite-size shift of T c behave 
like L~ l l u . The behaviour AT C (L) oc \jL d l 2 that would be naively dominant is destroyed 
by the disorder. The case of weakly bond-disordered spin models above the upper critical 



dimension needs a very careful analysis, as shown in [161 ]. 

To the best of our knowledge, the distribution of the pseudo-critical temperature in finite size 
spin glass models has not been studied numerically before: this is the object of the present 
note. Recent analytical work has predicted AT C (N) oc 1/iV 2 / 3 (where N is the number of 
spins, i.e. the system volume) for the (mean field) Sherrington-Kirkpatrick model (SK) for 
spin glasses [17|: we establish in this note that the realised scenario is indeed different. Very 
recently, while this work was being completed, ref. [3] has also tried (and failed) to verify 
numerically the analytical predictions of 17|. A former attempt to analyse numerically the 



distribution of the pseudo-critical temperature in the SK model was useful to investigate 
the numerical techniques of choice Ji|. 

Here, we present numerical results both for the three-dimensional (3d) and the mean-field SK 
spin glass models. In the 3d case, we show that the probability distribution of pseudo-critical 
temperatures verifies finite size scaling. From the scaling of this distribution we obtain a 
precise estimate of the critical temperature and of the critical exponent for the correlation 
length, v. On the other hand, we find that for the mean- field spin glass AT C (N) oc 1/iV 1 / 3 (in 
agreement with analytical findings for the scaling with N of disordered-averaged quantities 
in mean field models 20]). Since this is in plain contradiction with the results of Ref. 



we briefly revisit their analytical argument and show where the error in 17| stems from. We 
also believe that a second analytic conclusion of 13], stating that the AT C (N) is distributed 
according to a Tracy Widom probability law, is based on very shaky grounds, and we will 
give hints of the fact that it is not substantiated numerically. 

Our first step is to define a pseudo-critical temperature for a given finite-size sample. Ran- 
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dom bond and site diluted models allow 
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13) a straightforward definition of T C J (L) as the 



location of the maximum of the relevant susceptibility %■ 111 our case the situation is more 
complex (even if, as we will see, the analysis of the spin glass susceptibility xsg will be very 
useful and revealing). Here the relevant diverging quantity is 21] xsg 

XSG = J2^ S vh ~ ( S *)j(Sy)jf , (1) 
x,y 

(where the S x are the local spin variables) that is of order N in the whole low temperature 
phase, xsg is a continuously decreasing function of the temperature, and has no peak close 
to T c : this requires, as we will discuss in the following, a slightly more sophisticated analysis 
in order to extract a pseudo critical temperature. 

An alternative and simpler procedure is very straightforward: let us introduce it first. We 
first assume (as done in |q 13 1 ) that for a given disorder sample J the finite-size scaling 
of an observable P of dimension ( is (P)j — F((T — T c J (L))L 1 / i/ ), where F(-) is a L 
and J independent finite-size scaling function: the whole disorder sample dependence is 
encoded inside the pseudo-critical temperature T/(L). This is in fact an approximation 



since the scaling function has a residual J dependence We next build dimensionless 

combinations of operators: we call them J (T, L), and they are build in such a way to scale 

as 

J (T, L) ~ G((T - T C J (L))L 1/U ) , (2) 

where G(-) is a T, L and J independent finite size scaling function. A familiar looking 
combination is the single sample pseudo-Binder cumulant B J = (q 4 ) j / {q 2 )"j (notice that 
this is defined for a given disorder realisation): the genuine Binder cumulant is defined as 
B = E((q A )) / E({q 2 )) 2 . For sensible choice of y the solution Ty of the equation J (Ty , L) = 
y, with a disorder independent constant y, is a proxy of the pseudo-critical temperature, 
namely Ty = T C J (L) + C y L~ x l v = T c + DyL" 1 ^ with an L and J independent constant C y . 
For example for a function O that in the infinite volume limit is zero in one phase and one in 
the other phase any constant y in the interval [0, 1] will do: it is wise, in order to minimise 
the corrections to scaling to choose a legitimate value for y such that Ty is typically inside 
the critical region (the value of Ty depends on y and on the choice made of a dimensionless 
combination J ). 

We are also able to use xsg f° r determining T/(L): in this way we are able to monitor a 
quantity that diverges in the infinite volume limit, and to use it to extract a pseudo-critical 
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temperature. The approach used to define T C J (L) in this case is based on the same technique: 
we compare the single sample spin glass susceptibility to a value close to the average spin 
glass susceptibility at the critical temperature on a given lattice size. This measurement is 
a good proxy for the direct measurement of the position of an emerging divergence. 
We have applied these ideas to the Edwards- Anderson model in 3d and to the SK model. 
We used existing data obtained by the Janus collaboration on systems with L = 8 to L = 32 
for the 3d EA model [4], and from 22] for the SK model with iV ranging from 64 up to 4096. 
In both cases the quenched random couplings can take the two values ±1 with probability 
one half. 

The layout of the rest of this work is as follows. In Sect. [Til we discuss our numerical 
methods, and we present our results for the Edwards- Anderson model. An analogous analysis 
for the SK model is presented in Sect. IIHI This study is complemented in Sect. IIII Bl 
with our analysis of the analytically predicted scaling for the distribution of pseudo-critical 
temperatures. We also present in Sect. IIII CI an analysis of the distribution function of the 
pseudo critical points. Finally, we give our conclusions in Sect. [IVl 



II. THE EDWARDS-ANDERSON 3d MODEL 

The Hamiltonian of the model is 

H^d = — ^] S x J x y S y , (3) 

x,y 

where the sum runs over the couples of first neighbouring sites of a 3d simple cubic lattice 
with periodic boundary conditions. The Ising S x spin variables can take the two values ±1 
and the couplings are quenched binary variables that can take the value ±1 with probability 
one half. 

In order to analyse the single sample pseudo-critical temperatures of Edwards- Anderson 3d 
systems we need to construct several dimensionless quantities. We define the Fourier Trans- 
form of the replica-field q x = S X S X (S a and S b are two real replicas, i.e. two independent 
copies of the system evolving under the same couplings, but with different thermal noise) 

m = J2 q*j kx , (4) 
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which we use to construct the two-point propagator 

G J (k) = (0(fe)0(-fe)) J , (5) 

where (■■ - )j denotes the thermal average for the sample J. Since the smallest momentum 
compatible with the periodic boundary conditions is \k\ = 2-k/L we define 

k[ 1] = [2-k/L, 0, 0) , kf ] = (0, 2vr/L, 0) , kf ] = (0, 0, 2ir/L) , (6) 

and 

G J (k 1 )= 1 -J2 GJ ( k ?)- ( 7 ) 

i 

Similarly, the second smallest momentum is given by (2ir/L, 2it/ L, 0) and by the two other 
possibilities: we use it to define G(k 2 ). 
We consider the following dimensionless quantities: 

1 



(8) 



2Lsin(7r/L) [G^^ 
B = F y J , 0) 



Bg = ' m 

B J ~ GJ{kl) (11) 
In Eq. [9] we have used the global spin-overlap, which is defined as 

(12) 

where N = L d is the total number of spins. Let us start by considering the sample-averaged 
observables (which we denote by dropping the super-index J). Up to scaling corrections, 
they do not depend L at the critical point, 

0(T c ,L) = y c + 0(L- a ) , with a > . (13) 

For each value of L we can search for the temperature T^ y such that 

0(T£ u ,L) = y. (14) 

Then, provided we are not very far from the scaling region (so that y is not too different 
from y c ), we expect that 

T c L y ^T c + A y L-^(l + B y L-») , (15) 
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FIG. 1. Computation of the pseudo-critical temperature using the Binder ratio defined in eq. 
Here y = 1.51, that is a good approximation to B(T C ), and we plot B J (T) for three samples with 
L = 32. The red curve has exactly one solution (in our range) for B J (T) = y, which defines its 
Ty =1 51 : this is the normal behaviour of this observable. The green sample has three solutions, so 
we pick the largest one (green dot). The blue sample has B J (T) > y: in this case we take the 
lowest simulated temperature T ~ 0.7 as an upper bound for the pseudo-critical temperature (this 
ignorance will not affect our estimate for the median). Error bars are included in the plot and are 
very small. 

where we have included the first corrections to scaling. 

We can use this same approach to define a single- sample critical temperature T C J . Let us 
choose a fixed value of y close enough to the value y c defined from the sample average at 
the critical point. For each sample we use cubic splines to determine Ty ,J such that 

J (T y L ' J ,L) = y . (16) 

For some samples the J (T, L) turn out not to be monotonic: there can be several solutions 
to this equation. In those cases we simply pick the largest solution. This process is illustrated 
in Figure [TJ 

The motivation for this choice is simple. The physical meaning of a pseudo-critical tem- 
perature is a characteristic temperature that separates the paramagnetic phase from the 
low-temperature one. Indeed, any temperature Ty , solving the equation J (T y , L) = y, 
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is a temperature where non-paramagnetic behaviour has already arisen. Therefore, only 
the largest Ty makes sense as a divider among both phases. In fact, the non-monotonic 
behaviour of J {T,L) may be due to other reasons (in particular, temperature chaos), un- 
related to the paramagnetic/spin-glass phase transition. In any case, our definition will be 
justified a posteriori, on the view of the simplicity of the emerging physical picture. 
The values of Ty have a very wide probability distribution. For a few disorder samples the 
solution of eq. ( fl6l) falls out of our simulated range of temperatures and we only obtain an 
upper or (less frequently) a lower bound (see the blue curve in Figured]). In this situation the 
arithmetic average of the Ty is not well defined: we consider instead the median temperature, 
that we denote by Ty . Since, by definition, the median does not change as long as the 
proportion of samples without a solution is less than 50%, this is a robust estimator in 
these circumstances (we are well below this limit for all the cases considered, the typical 
proportion being about ~ 1%). From fig. [T| it is also clear that the statistical uncertainty 
over the determination of Ty in a given sample is very small as compared to the size of 
sample to sample fluctuations. 

In analogy with the sample averaged case (TT5l) . we make the ansatz 

Ty{L) —T c + AyL" 1 ^ , (17) 

where we have ignored sub-leading corrections. A fit to this equation would, in principle, 
yield the values of T c and v. However, for a fixed value y, we do not have enough degrees of 
freedom to determine simultaneously v and T c . 

Following the approach of [5], we get around this problem by considering n values of y 
at the same time: this allows us to fit at the same time all the resulting Ty(L), with fit 
parameters {T c , v, A yi , . . . , A Vn }: in other words we force the Ty obtained for different y 
values to extrapolate to the same T c with the same exponent. This procedure may seem 
dangerous, since we are extracting several transition temperatures from each of the J (T, L), 
that are correlated variables. However, the effect of these correlations can be controlled by 
considering the complete covariance matrix of the data. 

The set of points {T^(L a )} are labelled by their L and their y: we have data for C = 5 
different values of L, with L\ = 8, L 2 = 12, L 3 = 16, L 4 = 24, L 5 = 32 (4000 samples in 
all cases but for L = 32, where we have 1000). We also select n values of ?/j in the critical 
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region . . ., y n ). The appropriate chi-square estimator is 

n C 

x 2 = EE - T < - A/A^l^ Km - T c - A y]L -^] , (is) 

i,i=l a, 6=1 

where o-( ia uM is the covariance matrix of the set of Ty, which we compute by a bootstrap 



approach 23j (each point of the set is identified by L and y): it is a block-diagonal matrix, 
since the data for different L values are uncorrelated. 

Thus far we have considered just a single observable O, for different values of y. Since the 
fitting function ((TTj) is the same for the different O we have selected, with common T c and 
v, and only amplitudes differ, we can include in the same fit data for the four dimensionless 
quantities (IBT fTTT) . considering several values of y for each. In order to simplify the notation, 
from here on we shall denote our set of points in the fit as {T^(L)}, taking a as labelling 
both the observable O and the height y (so that it will range from 1 to An). 
We can use the usual disorder averaged spin-glass susceptibility xsg to arrive at yet another 
determination of the single-sample critical temperature, with the definition 

XsgC#) = Xsg(T c ) y , (19) 

with y close to one, and we expect T£ to follow the same scaling behaviour of (ITT]) : we 
include the values of T£ in the global fit to the individual pseudo critical temperatures. In 
this case the pathologies that affect the analysis of other single sample observables are far 
less frequent. 

We show the results of this combined fitting procedure in Figure El where we have included 
data for the four dimensionless ratios £/L, B, Bq and R12, using three values of y for each 
ratio (we plot the three fits for the same observable with the same colour). We have also 
used the data for T^, with y = 1 (we select the value of T c reported in [24J). We have 
discarded the L = 8 data, which showed strong corrections to the leading scaling of ffTTj) . 
The best fit gives 

T c = 1.104(6) , v = 2.26(13) , (20) 

with x 2 — 32.6 for 37 degrees of freedom (giving a P value of 68%). This results nicely agree 



with the determination of 24j . 



T c = 1.109(10), v = 2.45(15). (21) 
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FIG. 2. Global fit for the median single-sample critical temperature computed with different 
observables. The parameters in the fit are the common extrapolation point T c , the common 
exponent v and an individual amplitude for each set of data points. We include in the fit the data 
for L > 12. 



Ref. [24j includes corrections to scaling, as in ( TT5l) . with uj = 1.0(1). Ref. [25] gives a 
comprehensive list of estimates for T c and v. In order to take these scaling corrections into 
account and to include the L = 8 data in the fit, we can therefore rewrite (TT7I) as 



f J a {L) ~ T c + A a lT x l v (\ + B a L~ 



(22) 



where again we use the same u parameter for all the observables and all values of y. Unfor- 
tunately, our numerical data are not precise enough to allow a reliable determination of u, 
v and T c at the same time (the resulting error in u would be greater than 100%). We have 
been able to check consistency of our approach by taking the values of v and uj from [24] and 
fitting only for T c and for the amplitudes, including now the data for L = 8. The resulting 
best fit gives T c = 1.105(8), with x 2 — 41.9 for 38 degrees of freedom (P value: 35%): this 
is a satisfactory check of consistency. 

The results we have discussed make us confident of the fact that our determination of the 
single-sample critical temperatures yields reasonable results. We can now take the analysis 
one step further and consider the width of the distribution of T C J . We consider the two 
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FIG. 3. Global fit to AT^, computed using £ J /L, B J and x'sg- 



temperatures T+ and T a such that 



P{T J a > T+) = 0.16 



P{T J a < T~) = 0.16 . 



(23) 



The value 0.16 is such that the temperature interval [T~,T+] defines the same probabil- 
ity as an interval of two standard deviations around the mean for a Gaussian probability 
distribution. We define the width AT^ as 

T+ - T 7 



(24) 



The simplest ansatz for the scaling behaviour of AT^ is 



AT£ c± A a L~ 1 /"(l + B a L~ u ) . 



(25) 



In principle we could repeat the global fitting procedure that we have applied to the medians 
T£. Unfortunately, not all the observables that we considered in Figure [2] can be used 
to analyse AT^, since the distribution of some of them is too wide, so that the critical 
temperature of too many samples falls out of our simulated range of T, and the width 
defined in eq. [25] is undefined. Because of that we analyse AT^ by only using the Ty derived 
from £ J /L, B J and Xsg- The corrections to scaling are now stronger than for the median 
Ty, so that we cannot obtain a good fit to leading order even if we discard the data for 
L — 8. Using once again as an input the critical exponents from 24J and fitting for the 
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amplitudes we obtain a very good fit with x 2 — 19.1 for 21 degrees of freedom (P value: 
58%). The results of the best fit are plotted in Figure [3) According to the ansatz of eq. [2, 
the width of Ty must be equal to the width of T C J (L) (since C y is J independent): the width 
of Ty should accordingly be y independent. Figure [3] shows that this is indeed the case for 
AT^ L , but less so for AT^: as mentioned before the disorder independence of the scaling 
function in eq. [2] is only approximate 11]. 



III. THE SHERRINGTON-KIRKPATRICK MEAN FIELD THEORY 

The Hamiltonian of the Sherrington-Kirkpatrick mean field model is 

where the sum runs over all couples of spins of the system, the Ising Si spin variables can 
take the two values ±1 and the couplings are quenched binary variables that can take the 
value ±1 with probability one half. 



A. The pseudo-critical temperatures 

Our analysis of the mean field SK model is very similar to the one we have discussed in the 
case of the Edwards- Anderson model. Here we have seven values of the system size N (of 
the form N = 2 P , with p ranging from 6 to 12: in all cases we have 1024 disorder samples 
but for N = 4096 where we have 256 and for N = 128 where we have 8192), that makes the 
fitting procedure stable. It is also of use the fact that in this case the value of the (infinite 
volume) critical temperature, T c = 1 is known exactly. 
We consider three definitions of T c : 

• the one based on the single sample Binder cumulant of eq. [91 B J {Tg) = B(T C ); 

• one based on the low order cumulant D = E((q 2 )) / E((\q\)) 2 and the single sample 
quantity D J = (q 2 )j/(\q\)l i.e. D J (T J D ) = D(T C ); 

• one based on the spin-glass susceptibility Xsg^y) = Xsg(T c ). 
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For a given value of N, we use for the left hand sides B(T C ), D(T C ) and Xsg{T c ) the val- 
ues measured for this same size, at T c = 1. We solve eq. ffl6l) by using a simple linear 
interpolation. 

Like for the EA model, in some cases eq. (ITBT) can have more than one solution. We again 
choose the largest solution, which in the case of SK turns out to be always the one closer 
to the infinite volume value T c . In a few cases, for small values of N, the equation has only 
solutions outside the range of temperatures that was used in the parallel tempering Monte 
Carlo simulation (0.4 < T < 1.30). We fix this problem again by basing our statistical 
analysis on the median of the distribution and on the definition of the width given by eq. EH 
It turns out that these pathological cases are less numerous for the SK model than for the 
EA model, and that the width given by eq. EH is always defined. 

In terms of the number of sites N of the SK fully connected lattice the ansatz of eq. [T71 
becomes 

f/(L) ~ T c + A y N- l/{ud ^ ] =T C + A y N~ 1/3 . (27) 

We show in fig. [Hthe data for the median of the distribution of T^ 7 as a function of 1/iV 1 / 3 , 
together with the results of two best fits. We first notice that the data are well compatible 
with the fact that in the N — > oo limit T c = 1. The first fit is a linear fit to the form 
= 1 + a/N 1 / 3 (with N > 256 and x 2 = 4.44 with 4 degrees of freedom). This is a good 
fit for the large systems, but it fails below N = 256. We also show a (very good) best fit 
including the next to leading corrections, with an exponent 2/3 (including all N values, 
X 2 = 2.07 with 5 degrees of freedom). The analysis of the data for Tg and T5 leads to 
the same conclusions: here however the leading term (oc 1/N 1 ^ 3 ) has a small coefficient and 
the effect of the next to leading term is stronger. In conclusion our data are in excellent 
agreement with an asymptotic 1 + 0(1/N lj/3 ) behaviour for the median of the distribution. 
We show in fig. [5] the width of the distribution of as a function of 1/iV 1 / 3 , together with 
the results of two fits, namely ATj = ci/iV 1//3 , and AT/ 7 = C\N~ 1 / 3 + C2iV~ 2 / 3 respectively. 
The data are well compatible with AT = in the limit iV — > 00, as expected. The leading 
order fit, including iV > 256, has a x 2 — 5.36 with 4 degrees of freedom. The two-parameter 
fit gives an excellent representation of the data (with a x 2 — 2.475 with 5 degrees of freedom) 
including now the iV = 64 and iV = 128 points. Very similar results are obtained for ATg 
and AT^. 

In conclusion our finite size, numerical analysis of the SK model strongly support an asymp- 



14 



0.99 



?^ 0.98 



0.97 



0.96 



Leading order 
Corrections 



0.1 



N 



-1/3 



0.2 



0.3 



FIG. 4. The median of the distribution of T£ as function of 1/iV 1 / 3 , together with the results of 
two best fits: a leading order fit, and a fit including a higher order correction (see text). 
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FIG. 5. The width of the distribution of as function of N 1 > 3 , together with the results of two 
fits: a leading order fit, and a fit including a higher order correction (see text). 



totic 0(N 1 / 3 ) scaling behaviour for the width of the distribution of the pseudo-critical 
temperatures. 
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B. Scaling with the system size and the stability of TAP states 



Our results are in contradiction with the claim made in [17J that the width of the finite-size 
fluctuations due to quenched disorder of the critical temperature of the SK spin glass scales 
like A _2//3 . We show here that we can give support to our numerical finding by means of a 
very simple scaling argument. 



The SG susceptibility can be computed from the TAP free energy 26[ as 

XSG = ^TrAT 2 , (28) 

where Ai is the Hessian of the TAP free energy at the relevant minimum. If we seat deep in 
the paramagnetic phase, the only relevant minimum of the TAP free energy is rrii = for all 
i — 1, 2, . . . , N. Note that in the pseudo-critical region, where (/3-/3 c )A 1 / 3 ~ 1, is not at all 
obvious that the such TAP solution is the relevant one (for instance a sub-extensive set of 
sites of size N a , with a < 1, could have non- vanishing m;, or maybe one could have for all 
sites | m j | ~ N~ a \ with a' > 0): the following discussion is relevant only in the paramagnetic 
phase and for system sizes so large that (/3 C — (3) ^> A -1 / 3 . 



It was shown some time ago 27J that, at = /3 C , the smallest eigenvalue of the Hessian at 
the fully-paramagnetic TAP solution is of order A~ 2 / 3 . In Ref. 17| it has been argued that 



W = (1 - Pf + -j^j-^P , (29) 

where is a random variable that, in the limit of large N, converges in distribution to 
a Tracy- Widom random- variable 28|. In particular, note that at the critical temperature 
(3 C = 1 the A~ 2 / 3 scaling is recovered. 

Now, only for the purpose of discussing the crudest features of the scaling laws, let us assume 
that xsg is dominated by the contribution of the smallest eigenvalue: 

Xs J g ~ tttJ— , (30) 

iVA min,J 

(31) 



A 1 / 3 2- ( 32 ) 



The analysis of Ref. 



17] is based entirely on eq. (j3T 
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Now, note that eq. (I32I) implies that interesting behaviour appears only when (1 — 0) < 
l/N 1 / 3 : this makes sensible to replace /3^^ with . At this point, an implication 
emerges for the scaling with iV of the average susceptibility. We have that 

X sg(T, AO = N^GiN^il - /?)) , (33) 

where the scaling function G has the form 

/+oo 1 
„ d$PTw($) (^T^F- (34) 

In the above expression Ptw( < ^ ) ) is the Tracy- Widom probability density function. 
Eq. (I34p is not acceptable for two main reasons: 

• the function G(x) in eq. (1341) is ill-defined, as the integrand has a non-integrable sin- 
gularity at $ = —x 2 ] 

• if one devises some regularisation procedure, dimensional analysis would indicate that 
G{x) ~ x~ A in the limit of large x. However, in order to recover the correct critical 
divergence Xsg°° ~ V(l ~~ &)■> one obviously needs G{x) ~ 1/x. 

The solution of these two caveats is of course in the fact that the initial assumption, xic ~ 
l/(A^A^ in j), is incorrect. The contribution of ~ N eigenvalues is crucial in order to recover 
the correct scaling behaviour G(x) ~ 1/x. Thus, the only lesson that we may take from this 
oversimplified analysis is that the J-dependent SG susceptibility will probably scale as 

X J SG = N^{[N^(1-P)}+^) . (35) 

This is exactly the scaling ansatz we made at the beginning, where ^ ^ is some random- 
variable that (in distribution) remains of order 1 in the large- N limit. This result is consistent 
with our numerical findings. 

Merely rewriting eq. fl3T|) as eq. fl32l) suffices to make it obvious that the asymptotic statement 
in Ref. is incorrect: the width of the distribution of the pseudo-critical temperatures 
scales with iV -1//3 . Indeed, if as done in Ref. [l^], one simply derives in eq. ( 13T]) with respect 
to (3 in order to get the maximum of the susceptibility, one finds that 1 — f3 c j ~ iV -2 / 3 . But 
at such value of /3 Ci j we have that (1 — f3) 2 ~ iV -4 / 3 jMj^^j, In other words, at the scale 
of (1 — f3) ~ 1/iV 2 / 3 eq. (|3TT) predicts an essentially constant behaviour, hence the supposed 
maximum of the susceptibility (recall that at such value of the fully paramagnetic TAP 
minimum is probably no longer the relevant one) has no physical meaning. 



17 



0.5 



0.4 



0.3 



0.2 



0.1 








-4 



-2 





x 



FIG. 6. Probability density functions for the Gaussian distribution and for the Tracy- Widom 
distribution for the Gaussian Orthogonal Ensemble, both with zero average and unitary variance. 

C. The probability distribution of the pseudo-critical inverse temperatures 



We have discussed in some detail about the features of the pseudo-critical, sample dependent 
temperatures in the mean field SK theory, and we have determined their scaling properties. 
The next step, that we present here, is in the study of their probability distribution. Un- 
fortunately, there is not any clean analytical prediction for the shape that this distribution 
function should take in the large- N limit. The only proposal known to us was put for- 
ward in ref. [l^]: when properly scaled, the pseudo-critical temperatures should follow a 
Tracy- Widom (TW) distribution, in our case for the Gaussian Orthogonal Ensemble (GOE). 
Unfortunately, as we have explained in Sect. IIIIBl the reasoning leading to that prediction 
is flawed (although, on the long run, the prediction itself could be correct). 
Lacking an analytical guidance, we will simply check whether our numerical data can be 
described by either a TW distribution, or by the ubiquitous Gaussian distribution. Within 
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our limited statistics and system sizes, the two distributions turn out to be acceptable (the 
Gaussian hypothesis fits slightly better our data, but a Tracy- Widom hypothesis is certainly 
consistent as well). Given the preliminary nature of this study, we shall restrict ourselves 
to the simplest determination of pseudo-critical temperatures, the one coming from the 
spin-glass susceptibility. 

Let us start by noticing that the difference of a Gaussian distribution Fg{<P) and a TW 
distribution for the GOE Ftw{4>) is indeed very small. We show in figure [6] both a Gaussian 
and a TW distribution with zero average and variance equal to one: it is clear that they 
are very similar. Some numerical values can be of help. In the case of zero average and 
unitary variance a Gaussian has a fourth moment equal to 3, as opposed to 3.165 for a TW. 
The Gaussian is symmetric and has zero skewness, while the TW distribution has a small 
asymmetry, with a skewness equal to -0.29, The Gaussian has a kurtosis equal to zero, while 
a TW has a kurtosis equal to 0.165. We will use these numerical remark at the end of this 
section for sharpening the outcome of our quantitative analysis. 

It is clear that in this situation, where the two target distributions are very similar, one has to 
keep under very strong control finite size effects, that could completely mask the asymptotic 
behaviour. It is important to notice that the effects we are looking at characterise not only 
the tails but the bulk of the distribution. 

Let us give the basic elements of our approach. We consider the value of the pseudo-critical 
inverse temperatures computed from the spin glass susceptibility for different iV values. We 
try to verify if a relation of the form 



(where is a random variable with null expectation values and unit variance, Gaussian or 
TW like) can account for our numerical data and, if yes, to determine the scaling behaviour 
of 7at and of a^. We now define the variable (3 C = (f3 c — j)/a. Let us assume that for 
a system of size N we have Kn samples, and therefore Kn pseudo-critical temperatures 
(3 C (N, s), s — 1, 2, . . . , Kn- The empirical distribution function (EDF) is then 



where 9 is the Heaviside step function. Note that the parameters a(N) and j(N) in Eq. [36] 
are unknown a priori. They will be determined through a fitting procedure (see Eq. [38] 



(36) 




(37) 
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FIG. 7. (top): a(N), defined in equation [36| versus A r ~ 1 / 3 , and our best fit including the first 
scaling corrections, (middle): 'y(N), defined in equation 1361 versus A r_1 / 3 , and our best fit. Here 
the fitting function only includes the leading term, since this form already gives a good value for x 2 ■ 
(bottom): collapse of the EDF for different iV values and of the theoretical distribution function, 
as described in the text. Here we show the case of a Gaussian distribution. 
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FIG. 8. As in figure [3 but for a TW probability distribution. 

below). In particular, they typically take different numerical values if we adopt the TW 
hypothesis, or the Gaussian hypothesis. 

We define the distance among the EDF and the theoretical distribution as 

1 ^ N 2 

D = -r^J2\ F &(N,s))-H N c (N,s))} . (38) 
K N f— ' L J 
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In order to get an error estimate we repeat the procedure for 1000 bootstrap samples for 
each value of N. In fig. [7] we consider the theoretical hypothesis of a Gaussian distribution. 
We show in the top part of the figure a(N) versus N~ 1 ^ 3 , and our best fit including the 
first scaling corrections. In the middle frame we show 'y(N) versus iV~ 1//3 , and our best fit. 
Here the fitting function only includes the leading term, since this form already gives a good 
value for \ 2 ■ In the bottom frame we show the collapse of the EDF for different iV values 
and of the theoretical distribution function, as described in the text. In fig. [8] we show the 
same data for the hypothesis of a TW distribution. 

In short, even if the Gaussian hypothesis is slightly favoured over the TW one, this analysis 
does not allow us to decide clearly in one sense or in the other one. The values of D are 
always comparable among the two cases. The estimates of j(N) are indeed more consistent 
for the Gaussian case, but the difference of the quality of the two fits does not allow us to 
say a clear, final word. 

In order to try to sharpen our analysis we have used the Cramer- von Mises criterion 291-131]. 
We will not give here technical details (see, however, endnote 32]), but will only discuss 
the most important features and the results. When we normalise the distribution to zero 
average we and unitary width we introduce a correlations about the values to be tested: also 
because of that we find a better fit to our needs the two sample formulation of the criterion, 
with a non-parametric approach, where we have to start by fitting the test statistics (since 
we cannot use tabulated values, because we are determining a at and 7^ in Eq. [36] from 
our finite-size statistics). Again, as in our previous analysis, tests do not allow us to select 
a Gaussian or a TW distribution: they are both characterised by very similar levels of 
significance. 

A very simple analysis is maybe the most revealing. As we have discussed at the start of 
this section the fourth moment of a normalised Gaussian is equal to three, while the fourth 
moment of a TW distribution is equal to 3.165. We plot in fig.[9]the measured fourth moment 
of the probability distribution, both for data from xsg and for data from the pseudo Binder 
parameter, versus 1/N. The thick straight line is for a Gaussian distribution (where the 
value is three), while the thinner straight line is for a TW distribution. Again, these data 
do not allow for a precise statement, but they seem to favour the possibility of a Gaussian 
behaviour (the data for xsg gi ye maybe the clearer indication). 

Our conclusions is that, given the quality of our data and the sizes of our thermalised 
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FIG. 9. Measured fourth moment of the probability distribution, both for data from XSG an d 
for data from the pseudo Binder parameter, versus 1/N. Thick straight line is for a Gaussian 
distribution (where the value is three), while the thinner straight line is for a TW distribution. 



configurations, that do not go beyond N = 4096, a Gaussian distribution is favoured, but 
we cannot give a clear, unambiguous answer. 



IV. CONCLUSIONS 



We have presented a simple method to study the probability distribution of the pseudo- 
critical temperature for spin glasses. We have applied this method to the 3d EA Ising spin 
glass and to the fully connected SK models. Our results are in excellent agreement with a 
median of the distribution that behaves asymptotically like T c +0{L" l / u ) (or 1+0(N~ 1 ^ 3 ) for 
the SK model), and a width of the distribution that behaves like C(L" 1/iy ) (or 1 + C(A^" 1/3 ) 
for the SK model). The value of v we find for the EA model is compatible with state of 
the art results. Furthermore, even if our number of samples is modest as compared with 
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Ref. 



24| , our determination of v and T c is competitive. An analysis of the probability 



distribution of the pseudo-critical inverse temperatures for the SK mean field model does 
not lead to firm conclusions, but hints to a Gaussian behaviour. 
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